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1. Motivations and objectives 

The dynamic subgrid-scale (SGS) model (Germano et al ., 1991; Lilly, 1992) has 
proven successful in the large-eddy simulation (LES) of several simple turbulent 
flows, e.g., in homogeneous, incompressible flow with passive scalars and homo- 
geneous, compressible flow (Moin et a/., 1991); in transitional and steady plane- 
Pouiseille channel flow (Germano et a/., 1991); and in passive scalar transport in 
channel flow (Cabot, 1991; Cabot &; Moin, 1991). The dynamic SGS model, using 
eddy viscosity and diffusivity models as a basis, determines the spatially and tem- 
porally varying coefficients by effectively extrapolating the SGS stress and heat flux 
from the small, resolved scale structure, thus allowing the SGS model to adapt to 
temporally varying flow conditions and solid boundaries. In contrast, standard SGS 
models require tuning of model constants and ad hoc damping functions at walls. 

In order to apply the dynamic SGS model to more complicated turbulent flows that 
arise in geophysical and astrophysical situations, one needs to determine if the dy- 
namic SGS model can accurately model the effects of subgrid scales in flows with, 
e.g., thermal convection, compressibility, and rapid uniform or differential rotation. 

Theprimary goal of this work has been to assess the performance of the dynamic 
SGS model in the LES of channel flows in a variety of situations, viz., in temporal 
development of channel flow turned by a transverse pressure gradient and especially 
in buoyancy-driven turbulent flows such as Rayleigh- Benard and internally heated 
channel convection. For buoyancy-driven flows, there are additional buoyant terms 
that are possible in the base models, and one objective has been to determine if 
the dynamic SGS model results are sensitive to such terms. The ultimate goal is 
to determine the minimal base model needed in the dynamic SGS model to provide 
accurate results in flows with more complicated physical features. In addition, a 
program of direct numerical simulation (DNS) of fully compressible channel con- 
vection has been undertaken to determine stratification and compressibility effects. 
These simulations are intended to provide a comparative base for performing the 
LES of compressible (or highly stratified, pseudo-compressible) convection at high 
Reynolds number in the future. 

2. Accomplishments 

» 2.1 Large eddy simulation of time- dependent channel flow 

The dynamic SGS model was used in the LES of fully turbulent channel flow 
driven by a uniform streamwise ( x ) pressure gradient that is suddenly turned by 
a transverse (z) pressure gradient 10 times larger. The DNS of this case was per- 
formed by Moin et al (1990). They found, counterintuitively but consistent with 
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experimental results of three dimensional boundary layers, that the turbulence ki 
netic energy and shear production rate initially decrease and later recover. Until 
Durbin (1992, and in this volume), no Reynolds averaged type model had been able 
to reproduce this behavior. 

The LES was performed with a spectral- Chebyshev code (Kim et al . , 1987) on a 
32 x 65 x 32 mesh in a 47r x 2 x 47r/3 box (in units of channel half-width 6). The 
dynamic SGS model used a ratio of test to grid filter widths of 2 in the horizontal 
directions (using a sharp spectral-cutoff filter) and 1 in the normal ( y ) direction (i.e., 
no explicit filtering in y). Defining the effective filter width as A = (A* AyA*) 1 / 3 
gives a test to grid effective filter width ratio A/A = 2 2/3 . A Smagorinsky (1963) 
eddy viscosity base model was used whose coefficient, assumed to be a function of y 
and time, was calculated at each time step by averaging over horizontal planes (see 
Cabot, 1991). An ensemble of temporally developing flows was approximated by 
initially generating 15 fully developed turbulent channel flow fields separated in time 
by a sufficient amount to make them statistically independent. The initial channel 
flow fields were developed for a friction Reynolds number ( Re T — u TO S / u, where u TO 
is the initial friction speed and u is the molecular viscosity) of 180. The 15 fields 
were simultaneously advanced in time from t = 0 to 1.2 (in units of <5 /tt ro ), and 
statistics were generated for each field every A t of 0.15 and averaged together. The 
statistics from this LES were in good qualitative and quantitative agreement with 
those from the DNS (Moin et al., 1990), although the recovery in the turbulence 
kinetic energy in the LES occurred at a slightly later time than in the DNS. 

To test if it was the SGS model that was responsible for these good results in the 
LES or if it was due merely to an accurate portrayal of the large-scale interactions, 
a DNS was computed on the same coarse grid. The initial fields for the time- 
dependent calculation were first run to statistical equilibrium on the coarse grid, 
rather than simply turning off the SGS model in the LES initial fields, in order to 
avoid spurious transients due to the sudden drop in effective viscosity. The initial 
statistics for the coarse DNS and LES cases are thus not the same. The results of 
the coarse DNS were for the most part found to be in qualitative agreement with 
the well resolved DNS results of Moin et al. (1990), but the quantitative agreement 
was substantially poorer than was found using the dynamic SGS model. Thus much 
of the “three-dimensional” response of the turned channel flow is contained in the 
large-scale interactions, but the finer details require the SGS model. The greatest 
disagreement was found in the temporal behavior of the total (resolved and SGS) 
dissipation rate (Figure 1), which is to be expected since it depends to a larger 
extent on the different treatment of the small scales. In the DNS of Moin et al. 
(1990) and the LES, the dissipation rate has a complicated behavior near the wall, 
initially decreasing at the wall but increasing farther out in the near-wall region; 
the wall dissipation eventually begins to recover at t = 1.2. In the coarse DNS, 
however, the dissipation rate (which begins at a substantially higher level at the 
wall than in the LES) decreases both at the wall and in the near-wall region with 
no sign of recovery at t — 1.2. Such inaccuracies in the energy rates likely lead to 
the quantitative discrepancies in the velocity statistics. 
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FIGURE 1. Total dissipation rates near the wall (plotted as the distance from the 
wall in units of S) for channel flow turned by a transverse pressure gradient, (a) 
LES using the dynamic SGS model; (b) coarse DNS computed on the same grid. 

2.2 Large eddy simulation of thermal convection 
2.2.1 Base models 

Simple eddy viscosity and diffusivity SGS models, with some near-wall correc- 
tions, are commonly used in the LES of thermal convection (see Nieuwstadt, 1990, 
for a recent review). Some modelers employ additional buoyancy corrections (e.g., 
Eidson, 1985; Mason, 1989; Schumann, 1991). The eddy viscosity and diffusivity 
models that I have used to date a s the basis for the dynamic SGS procedure can be 
generalized in a form similar to Schumann’s (1991) “first-order” SGS model, which 
is a vast simplification of more general, second-order, Reynolds-stress-like equations. 
The model for the residual SGS Reynolds stress at an arbitrary filter level is 

t - ±Tt(t)I = — 2t/|S = -2C„A 2 aS , (1) 

where I is the identity tensor, C v is the coefficient of the eddy viscosity i/*, A is the 
effective filter width, and S is the strain rate tensor; a is a scale rate defined below. 
The residual heat (or scalar) flux is modeled by 

h = -C Q A 2 aB • VO , 5 = 1 + c 2 flV9/(a 2 - c 2 N 2 ) , (2) 

where C Q is the eddy diffusivity coefficient, 6 is the potential temperature, ft is the 
buoyancy vector (gravity times thermal expansion coefficient), and N 2 = fl * V0. 
The scale rate a is given, from SGS energy production = dissipation arguments, by 

a 2 = | [S 2 + (ci + c 2 )iV 2 ] + 1 1 [S 2 + (cj - c 2 )JV 2 ] 2 + Cl c 2 jV^ j 


(3) 
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where S 2 = 2S : S and = (/? * ^)(Vl9 • V0). The constant or coefficient c x is, 
in principle, C a /CV = 1/Pr*; c 2 is, in principle, related to the ratio of turbulent 
time scales of the velocity and potential temperature. Notice that (3) reduces to 
the normal Smagorinsky model scaling {a = S) for no buoyancy ( /? = 0) and 
that o’ 1 and a 2 — c 2 iV 2 are positive semi-definite if CiC 2 > 0. Also notice that the 
residual heat flux in (2) is anisotropic with respect to V0 for finite c 2 and /?, being 
enhanced in the direction of buoyancy forces. (Analogous anisotropic terms could 
be included in (1) by replacing S by B ■ S; but Schumann (1991) found that they led 
to realizability problems in his LES and so advocates dropping them.) For c 2 = 0, 
we can identify h with — where a* is the eddy diffusivity. 

The differences in the base models arise from different treatments of c x and c 2 ; 

A. The “scalar” model has c x = c 2 = 0. C u and C Q are determined as functions of y 
and t by the dynamic test-filtering procedure. This is the model for the dynamic 
SGS model employed by Moin et al . (1991) and Cabot & Moin (1991). I have 
applied it to Rayleigh-Benard convection. 

B. The “buoyancy” model has c x as a coefficient equated consistently with C Q IC U = 
1/Pr* and c 2 = 0 (isotropic eddy diffusivity). This requires an iterative solution 
of the eddy coefficients (Cabot, 1991) with a Newton’s (secant) method. It has 
been applied to Rayleigh-Benard convection and low-Pr internally heated channel 
convection. 

C. The “Eidson” model, after Eidson’s (1985) SGS model, is the same as B but 
with ci taken as a constant (2.5) corresponding to his best value of Pr< = 0.4 for 
the LES of Rayleigh-Benard convection. C u and C Q are determined, as in model 
A, with the dynamic procedure. I have applied this model to internally heated 
channel convection. 

D. The “Schumann” model has c x and c 2 taken as constants (2.5 and 3.0, respectively, 
which are near Schumann’s (1991) best values for the LES of planetary boundary 
layers). C v and C Q are determined, as in model A, with the dynamic procedure. 
This model has been applied to high-Pr internally heated channel convection. 

2.2.2 LES of Rayleigh-Benard convection 

Large eddy simulations of Rayleigh-Benard convection were performed with a 
spectral-finite difference code (Piomelli et al. , 1987) with the dynamic SGS model 
using base models A and B, the same filters as described in §2.1, and a mesh 
of 32 x 63 x 32. The molecular Prandtl number Pr was taken as 0.71 (air), and 
Rayleigh numbers Ra = 8\0AO\8 3 /va (where A0 is the wall-to-wall mean potential 
temperature difference) of 6.25 x 10 5 , 2.5 x 10 6 , and 1 x 10 7 were considered with 
horizontal-to-vertical aspect ratios of 5, 6, and 7, respectively. 

The buoyant (B) base model was found to give very similar results to the scalar 
(A) base model without buoyancy production terms. This probably happened be- 
cause the buoyancy term is generally less than, or at best comparable, to the strain 
term in (3) for this flow and because even with a different scaling the dynamic 
eddy viscosities and diffusivities tend to adjust to a similar level. The dynamic 
SGS model with the buoyant base model typically required only 2 or 3 iterations 
to determine the eddy coefficients consistently; this still doubled the computational 
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FIGURE 2. SGS eddy coefficients and Prandtl number from the LES of Rayleigh- 
Benard convection with Ra = lx 10 7 and Pr = 0.71 using the dynamic SGS model. 


cost of the SGS model and, considering the little difference it made to the results, is 
probably not warranted. Occasionally the iteration scheme failed to find solutions 
at some planes, perhaps indicating that no real solutions existed. The scheme gave 
up after 10 iterations; but converged solutions were always found a few time steps 
later a s flow conditions changed. 

The SGS eddy viscosity and diffusivity using base model B are shown in Figure 2 
with respect to their molecular values for the Ra = 1 x 10 7 case. Their fairly 
low values (of order 1 in the core) are a result of trying to resolve a reasonable 
amount of horizontal small scales near the wall. The dissipation due to the SGS 
model is comparable to that from the large scales in the core of the flow but becomes 
negligible near the wall. In fact, the eddy viscosity usually has small negative values 
in the viscous boundary layer though this has virtually no effect on the convective 
flow; it is not known if this is a real physical feature or an artifact of the poor 
horizontal resolution there. In contrast, the heat flux carried by the SGS model 
terms is negligible in the core of the flow but typically 20-30% of the total near the 
walls, which will affect the heat flux statistics. A concern is that the test filtering in 
the dynamic SGS model may not give accurate results near the wall since it usually 
samples in the energy-bearing part of the energy spectra there. The SGS Prandtl 
number is also shown in Figure 2. It is less than the standard value of about 0.4 
(Eidson, 1985) in the core, where I find values of 0.20-0.25, but it becomes larger 
near the walls, reaching 0.6. Sullivan &: Moeng (1992) found qualitatively similar 
results for Pr* in an a priori test of a DNS field but at levels 3-4 times higher. They 
used, however, an effective filter width ratio of 4 (versus my 2 2 / 3 ) and Pr = 1; they 
also revamped the dynamic procedure in a way that gives only positive values of 
so a direct comparison is difficult. 
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Large-scale statistics (such as rms velocity and potential temperature fluctuation 
intensities and velocity-temperature correlations) were found to be in good agree- 
ment with experimental measurements in air by Deardorff & Willis (1967) and 
Fitzjarrald (1976) and with previous LES results by Eidson (1985). The Nusselt 
numbers (Nu = 2£|V0| U ,/A0) of 7.7, 12.0, and 18.0 found for Ra = 6.25 x 10 5 , 
2.5 X 10 6 , and 1 x 10 7 using the scalar (A) base model are about 5-10% higher than 
the experimental values reported by Fitzjarrald (1976) (Nu « 0.13Ra° 30 in air) and 
Threlfall (1975) (Nu « 0.178Ra°' 280 in gaseous helium). A DNS for Ra = 6.25 x 10 5 
with the same code gave Nu = 7.2. A coarse DNS needs to be performed for one or 
more of these cases to determine the actual extent to which the SGS model improves 
the results. 

2.2.3 LES of internally heated channel convection 

Turbulent channel convection in water (Pr & 6) with uniform volumetric heat 
sources and cooled, no-slip walls has been examined experimentally by Ivulacki & 
Goldstein (1972) and numerically by Grotzbach (1982). This flow is asymmetric 
about the midchannel: the upper part of the channel is convectively unstable and 
the lower part is stable. The convective heat flux in the fully developed flow is 
typically downgradient in the exterior regions and countergradient in the interior. 
Because of this inherent asymmetry, the LES of this flow is expected to be more 
sensitive to the SGS model; it also allows us to test the behavior of the dynamic 
SGS model in transition from unstable to stable regions. 

Large eddy simulations were performed with a spectral- Chebyshev code (Kim 
et a/., 1987) for Pr = 0.2 at Ra = 1.25 x 10 5 on a 32 x 65 x 32 mesh and at 
Ra = 1.25 x 10 6 on a 32 x 129 x 32 mesh, and for Pr = 6.0 at Ra = 1.25 x 10 5 on 
a 32 x 65 x 32 mesh. Here Ra = \fi\ q6 5 /a 2 v, where q is the thermometric heating 
rate. All simulations used a horizontal-to- vertical aspect ratio of 4. 

For the low-Pr runs, I used both the scalar (A) and buoyant (B) base models in 
the dynamic SGS model. Although there were some differences in the v t and a t 
profiles for the low-Ra runs using different base models, the large-scale statistics 
were not particularly distinguishable. They shared the traits of having negative 
values of v% and/or a t near the walls; and Pr* had values of 0. 1-0.2 in the upper 
convective region, growing to values near unity near the unstable upper wall and the 
lower, stable region. Nusselt numbers at the upper wall were found to be about 5% 
greater than in DNS results (0. Hubickyj &; W. Cabot, unpublished). The profiles 
of Vt and with respect to molecular values and Pr* are shown for the high-Ra 
case in Figure 3 using the buoyant (B) base model in the SGS model. Except in 
the narrow viscous boundary layers, v t and a t are positive. In the core convective 
region (y/6 = -0.25 to 0.75), Pr* is about a constant 0.2 but grows to values of 1-2 
in the near-upper-wall region and the stable lower channel. The eddy diffusivity 
remains positive throughout the center of the channel where the large-scale heat 
flux is countergradient; this means that the SGS heat flux is downgradient in this 
region, counter to the large-scale flow, and that at acts rather to dissipate thermal 
fluctuations. Since the vertical temperature gradient is small in the central region, 
however, the SGS heat flux is negligible there and only becomes significant in the 
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FIGURE 3. SGS eddy coefficients and Prandtl number from the LES of internally 
heated channel convection with Ra = 1.25 x 10 6 and Pr = 0.2 using the dynamic 
SGS model with the buoyancy base model. Pr*. 

near-upper-wall region, attaining 20-30% of the total heat flux as in the LES of 
Rayleigh-Benard convection. The Nusselt numbers for this case are found to be 
about 10% higher than DNS results. (The large-scale statistics were again found to 
be fairly insensitive to the base model used.) 

The LES with the buoyancy base model experienced significant iteration prob- 
lems, most noticeable in the low-Pr, high-Ra run. Not only were there instances of 
failure to converge to a solution at some planes, more disturbingly there were clear 
instances when more than one solution existed and the values to which Pr* con- 
verged depended on the initial guess. (I needed to average the initial guesses over 
adjacent planes to get reasonable answers.) On the other hand, the LES with the 
scalar base model gave a broad drop in v t in the upper convective region, in poor 
agreement with the previous model (see Figure 4). Better agreement was found 
using the Eidson (C) base model, which includes the buoyancy production term in 
a less consistent but cheaper way than the buoyancy base model. The choppiness 
in v t in Figure 4 may be due in part to some numerical instability from advancing 
the SGS terms explicitly in the code at too large a time step, but it may also stem 
from filtering only in planes and not in the vertical direction, which would probably 
smooth the results considerably. 

For direct comparison with laboratory experiments, simulations with Pr = 6 
have been recently undertaken. The eddy diffusivities from the Ra = 1.25 x 10 5 run 
using the Eidson (C) base model are shown in Figure 5a. The eddy viscosity and 
Pr* are found to be negligible everywhere since the velocity in this case is almost 
completely resolved. However, near the upper wall I find Pr* « 5-7 (comparable 
to Pr). The eddy diffusivity in this case does have negative values in part of the 
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FIGURE 4. SGS eddy viscosity from the LES of internally heated channel convec- 
tion with Ra = 1.25 x 10 6 and Pr = 0.2 using the dynamic SGS model with base 
model A (scalar), B (buoyancy), and C (Eidson). 

central, countergradient region. Results using the Schumann (D) base model are 
also shown in Figure 5; a t in this case is defined as — h • V0/V0 • V0. Some minor 
differences in the central, countergradient region are noticeable. The residual heat 
flux resulting from these two models are shown in Figure 5b. It is seen that the 
Eidson base model only contributes to the heat flux in the upper convective region 
where the temperature gradient is appreciable while the Schumann base model 
contributes to the heat flux farther into the central region and gives comparatively 
more heat flux in the upper convective region due to the additional buoyancy term 
in Equation (2). Note that the SGS terms virtually vanish in the lower wall region 
where the flow becomes nearly laminar and that the dynamic SGS model allows 
a smooth transition between the turbulent and laminar regions. The LES results 
again tend to overestimate the Nusselt numbers by about 5% compared to DNS 
results; preliminary results indicate that coarse DNS computed on the same grid 
as LES overestimates Nu by more than twice as much. There is some discrepancy 
between experimental results (Kulacki & Goldstein, 1972) and numerical results 
(see Grotzbach, 1982), the former tending to give smaller Nusselt numbers and 
larger mean potential temperatures, the latter shown in Figure 6. The two different 
DNS results agree well but lie well below the experimental results; the LES results 
lie slightly below the DNS results (which make a fairer comparison). 

2.2.4 Conclusions from. LES results 

The dynamic SGS model has been used in the LES of a number of buoyancy- 
driven flows with different eddy viscosity /diffusivity base models that do or do not 
include buoyancy terms. I tentatively conclude from the results so far that the 
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FIGURE 5. SGS (a) eddy diffusivity and (b) vertical heat flux from the LES of 
internally heated channel convection with Ra = 1,25 x 10 5 and Pr = 6 using the 
dynamic SGS model with base model C (Eidson) and D (Schumann). 


buoyancy base model, which requires the consistent (iterative) determination of 
Pr*, is too computationally expensive and sometimes has either no real solution or 
multiple solutions. The “Eidson” base model, which simply sets Pr* to a constant in 
the model scaling, seems to provide a cheaper alternative that generally reproduces 
the buoyancy model better than the scalar model. It is not clear yet that the 
“Schumann” base model confers any real advantage over the others although it can 
accommodate, in principle, the countergradient heat flux that occurs in internally 
heated channel convection. 
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FIGURE 6. Mean potential temperature for internally heated channel convection 
with Ra = 1.25 x 10 5 and Pr = 6. ■ experimental data (Kulacki & Goldstein, 1972), 

o DNS (Grotzbach, 1982), DNS (0. Hubickyj & W. Cabot, unpublished), 

LES with base model C, and LES with base model D. 


2.3 DNS of fully compressible convection 

Direct numerical simulations of fully compressible, internally heated channel con- 
vection were performed using a fourth-order, explicit, finite-difference code (Thomp- 
son, 1990, 1992a,b). Simulations were performed for several different density and 
temperature stratifications at Ra = 1.23 x 10 5 (defined at midchannel) and Pr = 0.2 
in a linearly varying gravity. Fixed temperature, no-stress (free-slip) boundary con- 
ditions are used at the walls. The no-stress, impermeable walls are meant to ap- 
proximate free boundary conditions. For uniform volumetric heating rates, a mesh 
of 96 x 33 x 96 and horizontal- to- vertical aspect ratios of 4 or 5 are used; for uniform 
specific heating rates, a mesh of 64 x 65 x 64 and horizontal-to- vertical aspect ratios 
of 3 or 4 are used. 

The mean potential temperature profile from a low stratification, low Mach num- 
ber run was found to agree very well with the Boussinesq results of Cabot et al. 
(1990) for nearly the same values of Ra and Pr. For moderate to large density 
stratifications (central to wall ratios of a few to greater than 10) and moderate tem- 
perature stratification, the convection was found to be weaker due to the increase in 
viscosity and diffusivity (with inverse density) toward the walls; the Nusselt number 
was found to vary approximately as Nu — 1 oc (Ra 1 / 4 — Ra 1 / 4 )(/>„, /p c ) 3 / 4 , where 
Ra c ~ 1000 is the critical Rayleigh number for the onset of convection. The in- 
terior rms Mach number was found to be typically 0.20-0.25, increasing to about 
0.4 at the free-slip walls. Peak Mach numbers were found to be about 2.5 times 
the rms, and Mach numbers slightly in excess of unity were observed at the walls 
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in agreement with previous simulations by Malagoli et al. (1990). The compress- 
ible code did not require an additional high-order artificial damping built into it 
to compute these runs. Only weak shock features appeared to form because the 
high speed flows that form as the hot, rising interiors of convective cells expand 
horizontally along the walls tend to impinge on neighboring cells obliquely as the 
convergent flows plunge downward in cool, narrow downdrafts. Even a simulation 
with high temperature stratification (with central to wall ratio of ~ 30) with peak 
Mach numbers at the walls of 3.8 and occasional strong shock fronts was able to 
run a fair length of time without the artificial dissipation to damp two-delta waves, 
although it was eventually needed in this case. 

The levels of fluctuations in thermodynamic quantities relative to their mean 
values are found to be consistent with those of Chan & Sofia (1989) for simulations 
of deep stellar convection. As in their work, the rms pressure fluctuations were found 
to be almost equal to the turbulence kinetic energy everywhere in the convective 
region so that the relative pressure fluctuations scale as rms Mach number squared. 
An examination of the terms in the equation governing the potential energy P — 
p’ 2 /2~fp shows that they typically satisfy some of Zeman’s (1991) assumptions for 
a compressible boundary layer. The steady-state equation for P gives 


P 1 

— (u • VP + 27 PV * u + — u * Vp) — p'V ■ u' -p'u' * Vp 

p 7P 


Y JJt jyi 

-(p’u' * Vp' + 7 p ' 2 V - u') + (7 — 1) — — = 0 . 

7 P IP 


Here H is the net heating rate for the internal energy. As shown in Figure 7, term 
3 is a production term due to the pressure flux, which is very nearly balanced by 
the pressure dilatation in term 2. The remaining terms are higher order in Mach 
number squared and are negligible in moderate Mach number flows. Even in the 
high Mach number case cited previously, term 2 cancelled 60% of term 3. The pro- 
duction in term 3 is controlled here primarily by buoyancy terms since the pressure 
flux is proportional to the convective heat (enthalpy) flux and the pressure gradient 
is proportional to gravity from hydrostatic equilibrium. For convection then, un- 
like Zeman’s compressible boundary layer, the pressure flux should be modeled in 
terms of a thermal convection model, perhaps using the superadiabatic temperature 
gradient, rather than in terms of the normal density gradient. 

Compressional effects only appear to be significant at the (artificial) walls in 
the fully convective channels. Simulations with uniform specific heating rates are 
currently under way that feature convectively stable exterior regions bounding a 
convective interior. These should provide a better basis for determining compres- 
sional effects in the freely bounded convection; it is likely that acoustic effects are 
more important in the convectively stable exterior. We are also currently explor- 
ing whether the use of soundproofed, pseudo-compressible governing equations (like 
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Figure 7. Potential energy rates from equation (4) in fully compressible channel 
convection with high density stratification; pressure dilatation (term 2) and 

pressure flux production (term 3). 

Durran, 1989) in the simulations of highly stratified convection would be acceptably 
accurate and more efficient than the fully compressible simulations. 

3. Future plans 

3.1 “One- equation” local dynamic subgrid-scale models in channel flow 

Using locally defined coefficients from the dynamic SGS model has generally led 
to numerical instability due to persistent negative values of the SGS eddy viscosity. 
Ghosal, Lund & Moin in this volume (also see Wong, 1992) have proposed scaling 
the eddy viscosity with half of the trace of SGS residual stress ( k = r kk /2), which 
is evolved along with the flow. If the local k is driven to zero by negative eddy 
viscosities, the local eddy viscosity vanishes until k is replenished. This limits the 
duration of negative eddy viscosities and has been shown to stabilize calculations of 
homogeneous turbulence with local dynamic SGS modeling. We plan to implement 
this approach in channel flow. We also plan to implement Ghosal et al. s variational 
approach to determine the local dynamic coefficients consistently. 

An immediate problem arises in how to cast the fc-equation to give proper be- 
havior near and at the walls. Ghosal et al. use the form of a standard one-equation 
it- model with a SGS production term to evolve k at the grid-filter ( - ) level: 

Dk/Dt = u t S 2 + V • [(*/ + v D )Vk) - C E k 3 ' 2 / A , (5) 

where u t = CAk l/2 is determined by the dynamic test-filtering procedure (in which 
C is determined locally) and v D = CpAk 1 / 2 is the diffusive eddy viscosity. The 
constants or coefficients Co and C'e remain to be specified; they could be preset 
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constants or be themselves determined by a dynamic fitting procedure. First, k 
properly goes to zero at a no-slip wall as y 2 W (where y w is the distance from the wall). 
Since (5) is a second-order equation, only two boundary conditions are needed, 
namely k = 0 at either wall, but this generally results in k oc y w at the walls. 
Second, the uV 2 k term is generally finite at the wall, but it is difficult to make 
any other term in (5) balance it for plausible definitions of A, Cd, and Ce> (Note 
that u t in (5) always goes as y* w from the dynamic model.) Both of these problems 
can be addressed (but not necessarily solved) if we consider evolving the equation 
for q (where k = q 2 / 2) and understand the last term in (5) to be the model for 
the reduced dissipation rate e" = e - vVq • Vg, which goes as y 2 w at the wall. The 
additional term i/Vg • Vg must then be subtracted from vV 2 k in (5) to give ^gV 2 g, 
and the g-equation conveniently becomes 

Dq/Dt = cAS 2 + V * [{v + cjAg)Vg] + c^AVg • Vg - c e g 2 /A . (6) 

The lower case constants/coefficients in (6) differ from their upper case counterparts 
in (5) by various powers of y/2. Note that there is an additional source term in (6) 
from the diffusion term in (5). But now the term uV 2 q is generally finite and 
unbalanced at the walls (unless, e.g., c</A in the diffusive terms is made finite at 
the walls). However, even if this term is not balanced at the wall (in which case the 
numerical solution gives q — d 2 q/dy 2 = 0), one obtains the correct second-order 
asymptotic behavior for k at the wall ( k = dk/dy = 0). 

Tests of the sensitivity of the results to different treatments of the fc- or g-equation 
will need to be made for the LES of channel flow. A further modification that has 
been made to the channel flow code is to use top-hat (real space averaging) filters 
in the horizontal directions to assure that the trace of the residual SGS stress at 
the test-filter level is positive definite (which is not necessarily the case for spectral- 
cutoff filters). Also, top-hat filtering will also eventually be implemented in the 
direction normal to the walls for consistency with simulations of homogeneous flows, 
albeit not strictly commutative with normal derivatives on the stretched grid used. 
It also appears that a scheme must be developed for treating points with vanishing 
or negative q since equation (6) may have realizability problems. 

3.2 Further testing of the dynamic SGS model in thermal convection problems 

More large eddy simulations of Rayleigh-Benard and internally heated channel 
convection are needed to determine the optimal values of C\ and C 2 in equations 
(l)-(3) for the base models of the dynamic SGS model with plane averaging. The 
computational expense of computing them consistently at each time step with the 
dynamic test-filtering procedure is prohibitive (and ill-defined at some points), but 
sample calculations might be used to establish reasonable constant or functional 
values. For example, the value of Pr* = 1 fc\ is found to about a constant 0.2 in 
the core of several convective flows when C 2 = 0 (although this may be a function 
of filter sizes and the molecular Prandtl number). Some corresponding coarse-grid 
direct numerical simulations are also needed to gauge the effect of the SGS models. 
Filtering in the vertical direction, not explicitly done heretofore in the channel codes, 
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will also be implemented. Results from these volume-filtered channel simulations 
will be compared with previous DNS and plane-filtered LES results. The effect of 
including Leonard stress terms, similar to the mixed Smagorinsky-Bardina model 
(Piomelli et ai , 1987), will also be tested; these terms are generally non-dissipative 
but provide a fairly realistic level of local forward and backward scatter. 

More general base models for the subgrid scales are possible, especially in more 
complicated flows (e.g., with both buoyancy and rotation). Such models are being 
considered based on the governing equations for the residual Reynolds stress and 
heat flux, which closely resemble Reynolds stress equations for large-scale modeling 
(cf. Schumann, 1991). Dropping material derivative and diffusion terms, the gov- 
erning equations for residual stress (r), heat flux (h), and temperature intensity 
squared (k$) are 

At + t A* + 0h + h/? = n- 2 e, (7) 

r • V0 + A • h + /3ke = IT* — 2e$ , (8) 

h • V0 = -cm , (9) 

where A comprises the velocity gradient tensor and the mean rotation tensor (A,j = 
m } j — 2£Li€ijt) and A t is its transpose. The right-hand sides of (7)-(9) involve 
pressure-strain terms (FI) and dissipation terms (s) that must be modeled. A 
Smagorinsky model-like equation (1), for example, is recovered from (7) for stan- 
dard return-to-isotropy models of II and approximating the trace-free part of the 
left-hand side by 2r**S/3 • The importance of the more general terms will be tested 
in a priori tests of DNS data. The need for explicit rotational terms in the dynamic 
SGS base model will be tested with the LES of some rotating flows. A base model 
with rotational effects might be based on the above equations with rotation entering 
through A and/or through the models for II and e. 

The net amount of energy and dissipation that the dynamic SGS model can 
represent in channel flow has been limited due to the reduction of both normal 
and horizontal length scales near the no-slip walls, which causes the test filter to 
eliminate scales with a significant fraction of energy. Large eddy simulations for 
channel convection with no-stress walls will be performed in an attempt to improve 
on this situation. However, only flows with smaller scale disparity (freely bounded 
or with matching to near-wall solutions) will probably be able to use the dynamic 
SGS model efficiently at very high Reynolds numbers. Finally, we plan to implement 
the dynamic SGS model in our compressible convection simulations, starting with 
the form used in the LES of homogeneous compressible flow by Moin et ai (1991). 
The simulations now in progress with convectively stable exterior regions freely 
bounding the interior convective region should be more suitable for the dynamic 
SGS model by eliminating (or at least severely reducing) the amount of turbulence 
at the impermeable walls. 
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